% Define log-likelihood function for Hawkes process
function ll = hawkes_loglik(params, times)
    mu = params(1); % Background intensity
    alpha = params(2); % Excitation factor
    beta = params(3); % Decay rate
    
    N = length(times); % Number of events
    ll = -mu * times(end) + N * log(mu); % Background term
    
    % Excitation contributions
    for i = 2:N
        intensity = mu + alpha * sum(exp(-beta * (times(i) - times(1:i-1))));
        ll = ll + log(intensity);
    end
    ll = -ll; % Negative log-likelihood for minimization
end

% Initial parameter guesses
params0 = [0.1, 0.5, 1.0];

% Optimize parameters
options = optimset('Display', 'iter');
opt_params = fminsearch(@(params) hawkes_loglik(params, jump_times), params0, options);

% Display fitted parameters
disp('Fitted Hawkes Parameters:');
disp(opt_params);

% Extract parameters
mu = opt_params(1);
alpha = opt_params(2);
beta = opt_params(3);

% Compute intensity over time
time_grid = linspace(0, jump_times(end), 13200);
intensity = mu * ones(size(time_grid));

for i = 1:length(jump_times)
    intensity = intensity + alpha * exp(-beta * (time_grid - jump_times(i))) .* (time_grid >= jump_times(i));
end

% Plot
plot(time_grid, intensity, 'LineWidth', 1.5);
xlabel('Time (seconds)');
ylabel('Intensity');
title('Hawkes Process Intensity');
grid on;
